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ABSTRACT 

Astronomical optical interferometers (01) sample the Fourier trans¬ 
form of the intensity distribution of a source at the observation 
wavelength. Because of rapid atmospheric perturbations, the phases 
of the complex Fourier samples (visibilities) cannot be directly ex¬ 
ploited, and instead linear relationships between the phases are used 
(phase closures and differential phases). Consequently, specific 
image reconstruction methods have been devised in the last few 
decades. Modern polychromatic 01 instruments are now paving the 
way to multiwavelength imaging. This paper presents the derivation 
of a spatio-spectral (“3D”) image reconstruction algorithm called 
PAINTER (Polychromatic opticA1 INTErferometric Reconstruction 
software). The algorithm is able to solve large scale problems. 
It relies on an iterative process, which alternates estimation of 
polychromatic images and of complex visibilities. The complex 
visibilities are not only estimated from squared moduli and closure 
phases, but also from differential phases, which help to better con¬ 
strain the polychromatic reconstruction. Simulations on synthetic 
data illustrate the efficiency of the algorithm. 

Index Terms — ADMM, irregular sampling, phases estimation, 
proximal operator, optical interferometry 

1. INTRODUCTION 

The long-standing observation technique called Astronomical Inter¬ 
ferometry (AI) has lead to major discoveries during the last century. 
In optical wavelengths, the Very Large Telescope Interferometer in 
Chile will host in Summer 2015 two next generation instruments. In 
the radio domain, international efforts are being devoted to design 
and build a million-receptor array to be operational in 2020’s, the 
SKA (see https://www.skatelescope.org). 

In AI, the spatial position of each pair of receivers (telescopes or 
antennas) defines one of the N\> baselines of the telescope/antenna 
array. In absence of any perturbation and at a given observation 
wavelength, a pair of receivers with baseline b provides a complex 
visibility defined by y x = /(^-), which corresponds to a sample of 
the Fourier spectrum of the intensity distribution of interest (l). 

The sampling function is fully defined by the positions of the 
interfering receivers and by the observation wavelength. Despite the 
Earth rotation, which changes the geometric configuration of the ar¬ 
ray w.r.t. the line of sight and thus provides additional samples for 
observations acquired at different times, the sampling of the Fourier 
space remains extremely sparse in practice. AI is thus a typical in¬ 
stance of Compressed Sensing (CS), where the underdetermination 
of the reconstruction problem (restore an intensity distribution from 
few projections in the Fourier space) can be alleviated by exploiting 
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the intrinsic sparsity of astronomical scenes. This fact was early rec¬ 
ognized by Astronomers through the CLEAN algorithm (2j, which 
is a kind of Matching Pursuit assuming point sources (stars). 

In the last ten years, research in the field of interferometric 
image reconstruction (especially in radio) has been mainly driven 
by advances in sparse representation models and in optimization 
methods. Sparse representations using various redundant dictionar¬ 
ies have been used, including synthesis (3j and analysis 0 priors, 
and combination of both 0. 

This paper focuses in Optical Interferometry (01), for which the 
problem is more difficult than in radio. First, 01 does not allow ac¬ 
curate measurements of the phases (<p) of the visibilities because of 
rapid atmospheric perturbations. This excess of missing information 
(w.r.t. the radio case, where phases are available) can be partially 
recovered by means of linear relationships between the phases, the 
phase closures. This technique combines triplets of phases measured 
by different telescopes and produces a phase information which is 
theoretically independent of the atmospheric turbulence. Second, 
the number of telescopes in 01 arrays is much less than the number 
of antennas in radio arrays, whence a much sparser Fourier sampling 
in optical. 

In this framework, a classical and well-understood strategy for 
image reconstruction is to adopt an inverse problem approach, where 
missing information is mitigated, and hopefully compensated for, by 
a priori knowledge jlj. In this case, the image reconstruction algo¬ 
rithm aims at finding an intensity distribution that minimizes a cost 
function composed of a data fidelity term, which is related to the 
noise distribution, plus a regularization term and possibly other con¬ 
straints, which are related to prior knowledge. Following this path, 
various algorithms have blossomed in the last twenty years. Most 
of the proposed algorithms rely on gradient descent methods ( WIS - 
ARD (SJ, BSMEM (7), MiRA 0, BBM 0, IRBis (TO)). A different 
approach is used in MACIM JTlJ and in its evolution SQUEEZE GD’ 
which rely on Markov Chain Monte Carlo (MCMC) method. 

The case of polychromatic observations, which is under focus 
here, has recently been made possible through the advent of mul¬ 
tiwavelength interferometers. In this case, an astrophysical source 
is described by an intensity distribution which is a function of the 
wavelength, and the inverse problem aims at recovering the spatio- 
spectral (3D) distribution of the source. This objective adds to the 
intrinsic underdetermination of the 01 image reconstruction problem 
the computational complexity of solving a large scale inverse prob¬ 
lem. 

Polychromatic 01 image reconstruction has recently become a 
very active domain of research. A spatio-spectral image reconstruc¬ 
tion algorithm named PAINTER (for Polychromatic optic AI INTEr¬ 
ferometric Reconstruction software) with a publicly available source 
code has been proposed in ED- This algorithm relies on the visibili¬ 
ties (7 2 ) as measures of the source power spectrum and on two types 
of turbulence independent phases differences: the phase closures (ip) 



at each wavelength and the differential phases (14) , which are de¬ 
fined as the phases relatively to a reference channel. Those consti¬ 
tute additional turbulence-independent observables of the phases in 
multi wavelength observation modes. 

The objective of this communication is to present a number of 
improvements brought to the prototype version of PAINTER algo¬ 
rithm. These improvements i) regard the accuracy of the represen¬ 
tation model, which involves sparsity in analysis through union of 
bases and ii) make the algorithm highly scalable. The resulting 
changes represent an in-depth modification of the original version. 
A totally new version of the source code is also publicly available. 

The paper is organized as follows: section [2] introduces the no¬ 
tations and data model. In section[3]we present an detailed model of 
phase relationships. Section[4]tackles the inverse problem approach. 
Section [5] derives the resulting 3D reconstruction algorithm. Perfor¬ 
mances of the algorithm are presented and analyzed in section [ 6 ] 

2. DATA MODEL AND NOTATIONS 

Let y An be the complex visibility at the spatial frequency b/A n , 
and let y Xn be the column vector collecting the set of complex vis¬ 
ibilities corresponding to all available baselines at wavelength A n . 
The complex visibilities can then be related in matrix form to the 
parameters by the direct model EQD 

yXn = F A„ (1) 

where F Xn is obtained from a Non Uniform Discrete Fourier Trans¬ 
form (NuDFT) p~5) at the spatial frequencies imposed by the geom¬ 
etry of the telescope array and by the observation wavelength A n . 
Note that F Xn is not an orthogonal matrix. The previous expres¬ 
sion describes the complex visibilities by wavelength. A compact 
notation including all wavelengths and baselines is 

y = F x, F = ©^ F Xn (2) 


where F is a block diagonal matrix with each block referring to 
the NuDFT at a particular wavelength. Vector y concatenates the 
complex visibility vectors ( y Xn of Eq. [T]) for all wavelengths into a 
NbN\ x 1 visibility vector, with associated moduli 7 and phases <p 
given by 


y = [v XlT ,---, v X nxT ] , 7 = \v\, ^ v (3) 

In order to analyze the chromatic variation of the visibilities y Xn and 
of the images x Xn over the N\ wavelengths, we also introduce the 
Nb x N\ matrix Y and the x N\ matrix X defined as: 

X - vec -1 x , Y — vec -1 y (4) 

To clarify the use of a matrix notation note that the n th column of 
X, denoted as X n , corresponds to the vectorization of the image at 
the wavelength A n while the p th line describes the variation of pixel 
p along the wavelengths (i.e., a spectrum). 

3. MODEL FOR PHASE RELATIONSHIPS 

In the presence of atmospheric turbulence, the beams received at 
each telescope are affected by random and different optical paths, 
which corrupt the phases measurements of the complex visibilities. 


The “atmospheric corrupted” visibilities at a given wavelength A n 
for the base b a ,b involving telescopes a and b can be modeled as: 

y X a% = 7aH eX P (i [<Pa?b + Va" ~ Vb "]) ( 5 ) 

where p Xn is the uncorrupted phases and p Xn are perturbation terms 
related to telescopes a and b. To overcome the difficulty of phase es¬ 
timation, turbulence independent quantities need to be constructed. 


3.0.1. Phase closures 

The closure phase allows to get rid of atmospheric effects for triplets 
of complex visibilities. In presence of turbulent measurements, the 
closure phase (ip) is defined as the phase of the bispectrum (7j, i.e., 
the Fourier transform of the triple correlation. For three baselines 
b tt , b, bb ;C and b ac corresponding to a triplet (a, b : c ) of telescopes, 
the closure phase is defined as: 

_ / _ A n A n A n _ ,_A n 1 A n , _A n _ i X n , _ A n f 

t a,b,c - Z y a ,byb, c ya,c - Va,b + Pb,c “ ^a,c - ,b,c<P ( 6 ) 

where p Xn is the vector containing all unperturbed phases for wave¬ 
length A n , and h Xn b c is a sparse row vector with only three non zeros 
entries that take values {1,1, —1}. If N t denotes the number of tele¬ 
scopes, it is possible to show that (Nt — 1 )(N t — 2)/2 independent 
closure phases per wavelength are available (16). 


3.0.2. Differential phases 


For one baseline h a ,b, differential phases (Ap) measure the phase 
evolution in wavelength with respect to a reference phase channel. 
Because the phase turbulence terms on each telescope p Xk and p Xxti 
are, to a first approximation, independent of the wavelength ED and 
the differential phases defined by 


A Aef = z _ z y >w = ^ (7) 


I Aref — ,n^k 


A r ef _ l,A/;.,A re f 


are essentially not affected by the atmospheric perturbation. The 
reference channel can be chosen as one of the available channels. In 
this case, N\ — 1 independent differential phases are available per 
baseline. Without loss of generality, we chose Ai as the reference 
channel and Ap Xk b Xl = h Xk b Xl cp , where h Xk b Xl is a sparse row 
vector with only two non zeros entries that take values { 1 ,- 1 }. 


3.0.3. Model for all phase relationships 


The combination of the differential phases and phase closures into a 
global model will improve the phase estimation: indeed, the phase 
closures constrain the phases of a triplet of bases at a fixed wave¬ 
length, while the differential phases constrain the phases dependence 
in wavelength for a given base. To derive this model, we denote by 

= ( 1 ( N\ — 1) ® IjV b |I(AT^-J)xJV b ) 

the matrix concatenating all vectors h x ’ y Al of Eq.|7jin its rows. Sim¬ 
ilarly, = I Nx ® H ' is a block diagonal matrix that replicates 
the matrix H Xl concatenating the vectors h Xn b c of Eq. 6 in its rows. 
The information from the phase closures and differentiaTphases can 
then be collected in a global vector £: 
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- H i> 
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5 s — 
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( 8 ) 


where is the vector of all phase closures and A ip the vector of all 
differential phases. 








4. INVERSE PROBLEM APPROACH 

According to Eq.[ 8 ]and notations defined in Eq. [5] a data model for 
phases relationships and squared moduli can be written as: 

i = H(p + r) v £ = 7 2 + tJc (9) 

where 77 ^ and 77 ^ account for noise and modeling errors. Classical 
assumptions on their distributions are considered here. The noise 77 ^ 
is assumed to be jointly independent and Gaussian j 6 | and the noise 
77 ^ is assumed to be jointly independent and marginally Von Mises 
distributed 1 8 j|. Writing the opposite logarithm of the joint likelihood 
of £ and C leads to 

<? data (*) = a/(y 7 ) (10) 

where a and [3 are relative weighting terms and 


of X, which are the images at each wavelength processed indepen¬ 
dently. H s is a dictionary composed by the concatenation of the 
first eight orthonormal Daubechies wavelet bases (Dbl-Db 8 ) and a 
Haar wavelet basis. This type of regularization was recently used 
in radioastronomy |4j. H x X T operates on the rows of X to con¬ 
nect the pixels between wavelengths. In the present work H x im¬ 
plements a Discrete Cosine Transform (DCT) but this can be easily 
replaced by any union of orthogonal bases. Note that the related ma¬ 
trix H x is also an orthogonal matrix. Finally, g £ , fi s and g\ are 
hyper-parameters, which control the weights of the associated reg¬ 
ularization terms. Note that the previous version of PAINTER |13| 
was based on total variation regularizations. Besides the fact that this 
edge preserving prior does not perfectly match the smooth nature of 
the astrophysical sources, a major drawback of this choice comes 
from the heavy computational cost related to the non-orthogonality 
of the underlying image transform. 


5 C (y 7 ) = Hrr kn -7n) 2 > 7 = 1-^*1 (11) 

n 

9 i {y 4 ’) = COS (h m <P ~ £ m )> ¥> = ^ Fx (l 2 ) 

m 

Notations y 1 and are used above to underline that the first term 
depends only on the modulus and the second only on the phase of 
y, with y — Fx. The constant oj n is the variances of £ n . The 
constant K m is related to the variance of £ >rn by var(£ m ) = 1 — 
Ii(Km)/lo(Km) where Ij is the modified Bessel function of order 
j. For a given var(£ m ), K m is computed solving numerically this 
equation. 

As explained in the introduction, the problem is severely ill- 
conditionned owing to the poor coverage of the Fourier space. This 
requires tackling the image reconstruction as a regularized optimiza¬ 
tion problem (id We will adopt here an objective function of the 
form: 

x <— minimize (g data (x) + f eg (x)] (13) 

where the 3D image x can be constrained to have a spatially limited 
support II. Further constraints such as non negativity can be added 
in / reg (cc), which contains all the regularization terms. The support 
constraint is not included in f eg (x) for technical reasons related to 
the ADMM methodology described below. 

4.1. Regularizations and constraints 

01 images are by nature non negative and sometimes contain sources 
that are spatially localized. However, specifying the properties of the 
object parameters x only in terms of non negativity and spatial sup¬ 
port is usually not a sufficient prior. It follows that the use of regular¬ 
ization terms to emphasize some inherent a priori knowledge about 
the image structure is necessary. Following the matrix notation for 
the 3D object as defined in Eq. @ PAINTER in its current form in¬ 
cludes the ridge regularization, motivated by the poor conditioning 
of the NuDFT operator and spatio/spectral regularizations. The sup¬ 
port constraint is defined by the parameters space II in Eq. |T3] and 
the non-negativity constraint by the regularization term 1 R +(X). 
Consequently the regularization function in Eq.[T3| writes: 

r g (x) = 1 M+ (X) + ^-\\x\\l + fisWi^xiu + MA||Jr A x T || 1 

(14) 

H and H x are the matrices associated respectively with the spatial 
and spectral analysis regularizations Jl5) . H S X acts on the columns 


5. 3D RECONSTRUCTION ALGORITHM 


Owing to the unavoidable non convexity of the problem as defined 
by Eq. [l3](see e.g. in | 6 j), the vast majority of image reconstruc¬ 
tion algorithms rely on a descent optimization principle. So does 
PAINTER by using the flexibility of the Alternate Direction Meth¬ 
ods of Multipliers (ADMM) algorithm, which was already used in 
(B) to reconstruct stellar spectra of point sources from complex vis¬ 
ibilities. 

The optimization problem of Eq. 13 where g data (x) is given by 

Eqs. 


10 


12 


and the regularization term / reg ( x) is given by Eq. 


equivalent to: 


14 


minimize a g c {y 1 ) + /? g i (y <l> ) + ^-\\X\ 

y^,y<t>,y,x,Pen,T,S,V Z 


lf + 


s.t.: y 


y, y 


4> _ 


•••i m+ (P) + Ms ||t||i + M a||V||i 

y,y = Fx,T = H S X, V = SH X , S = X 


Auxiliary variables related to the complex visibilities: y, y 1 , 
y^ have proper Lagrange multipliers v y ,v^, v<f> and share the same 
augmented Lagrangian parameter p y . The auxiliary variables intro¬ 
duced by the regularization, P, T, V, S , have Lagrange multipliers 
Tp, Tt, Ty, T s and augmented Lagrangian parameters pp, pr 
and ps for V and S. The n th column of T. is denoted as v Xn . Min¬ 
imisation of the augmented Lagrangian leads to solve alternatively 
and iteratively the following steps: 

I. Minimization w.r.t. y 1 . Denoting y 1 = y + p y 1 v 1 

y 7+ = argmin ag c (y' t ) + \pyWy 1 - y 1 \\l 
y 1 Z 

This minimization is analytical. It comes down to find the real 
root of a cubic function using Cardano’s method. See (H) 

II. Minimization w.r.t. y (f> . Denoting y^ = y + py 1 v$ 

y* + = argmin /3/(j/) + \pyWy 4, - y^111 

This minimization is numerically solved using the limited mem¬ 
ory BFGS algorithm fl8) . 

III. Minimization w.r.t. y. 

y + = l(y 1+ + v^ + Fx- ppvy) 

This consensus step leads to complex visibilities reconstruction. 







IV. Minimization w.r.t. x. This step operates separately on each 
wavelength: 

C Xn X+ =F x " H (/ o y Y+ - v£") + H sT (p T T n - v x ") + 

• • • (p P P n - Vp"') + (ps s„ - 

With C Xn — p y F XnH F Xn H- 77 1, where 77 — p e -\-LpT + pp + 
ps . The Identity matrix in C Xn comes from the orthogonality 
of the L wavelet bases (L = 9 here) used for the spatial regular¬ 
ization. Computation of the right side term is realised using an 
inverse NuDFT and a inverse discrete wavelet transform. C Xn 
can be inverted using the matrix inversion lemma: 

c ' An_1 = ^(I - F XnH [F Xn F XnH + p-Vr^ An ) 

The number of rows in F Xri being small in 01 (small number of 
optical bases), the inner inverse matrix on the right side of the 
equality can be precomputed. Final computation of X+ is then 
realized using a NuDFT and an inverse NuDFT. 

V. Minimization w.r.t. P. Denoting P n = X+ + p^ > 1 2 3 4 v p Xn \ 

P+ = P R+ o P n (P) 
where Pc is the projection on the set C. 

VI. Minimization w.r.t. T. Denoting T n = H S X+ + p^}v t Xri '. 

T+ = Soft Ms (T n ) 

where Soft a (*) is the soft threshold operator. This step operates 
separately on each wavelength. 

VII. Minimization w.r.t. S. 

2 S+ = (V- Ps 1 r v )H x + X+ + p^Ys 

This step operates separately on each voxel. Note that in gd 
the right term of previous equality was multiplied by the inverse 
of H x H xt + I. The use of an orthogonal DCT reduces this 
term to 21 . 

VIII. Minimization w.r.t. V. Denoting V n — S+H X + Py 1 v v Xn : 
F+ = SofW(V) 

This step operates separately on each voxel. 

IX. The Lagrange multipliers are updated in the standard way, G3- 

A first acceleration of the proposed algorithm w.r.t. [131 relies on the 
use of orthogonal analysis dictionaries as detailed in the steps IV. 
and VII. The bottleneck of G3 was step IV. The computational time 
(matrix-vector multiplication in step IV.) is reduced from O (iV*) 
to O (4 N% log(N x ) + iVb(2 + iVb)) in the current version. 

Non-equidistant Fast Fourier Transform (NFFT, http://www. 
nfft.org) are used to efficiently compute the NuDFT and inverse 
NuDFT. Moreover, to take advantage of the possibility to parallelize 
the steps IV., VI. and VIII. w.r.t. the wavelengths or the voxels, the 
PAINTER algorithm has been totallv reimplemented in Juli^] This 
implementation relies on OptimPacljjfor the limited memory BFGS 
algorithm of step II., on the NFFT packag^] and on the discrete 
wavelet transform packag^] Optimisation in step II. is currently the 
bottleneck of the algorithm. The sources are publicly available at 
http ://w w w- n. oca.eu/aferrari/painter/ 

1 http://julialang.org 

2 https://github.com/emmt/OptimPack.jl 

3 https ://github. com/tknopp/NFFT.jl 

4 https ://github.com/JuliaDSP/Wavelets .jl 


6. SIMULATIONS 

This section presents simulation results obtained with realistic noisy 
synthetic data. The standard deviation (s.d.) of the noise for the 
phases is cr^ = 1 rd/SNR (in radians) and the s.d. of the noise for 
the amplitudes is cr 7 = 7 /SNR. SNR is set to 30 dB in both cases. 
The considered instrumental configuration is that of the 2004 Inter¬ 
national Beauty Contest in Optical Interferometry fl9) . The data are 
acquired in 13 acquisition epochs and in 8 equispaced wavelength 
channels in the range 1.47 qm — 1.56 pm. The resulting Fourier 
coverage, including the Earth rotation effect, is shown in Fig. [T] 
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Fig. 1. Spatial frequencies coverage 

The considered sources consist in two resolved stars for which 
the diameter and the brightness distribution vary in wavelength (see 
the left column of Fig.[2]). To initialize PAINTER we used for the 
N\ images the same image composed of a centered Dirac delta func¬ 
tion. The size of the reconstructed image at each wavelength is 
128 x 128 pixels.The algorithm was stopped after 1000 iterations. 

The estimated objects are shown in the right column of Fig [2] 
which shows that the shape and the diameter evolution of the sources 
are clearly well reconstructed. The relative mean square error of 
the reconstructed image is 7dB with G3 and 3.4dB with the pro¬ 
posed algorithm. The variation of a source’s integrated brightness 
as a function of wavelength is an interesting information per se and 
was thus also investigated (Fig. [3}. The integrated brightness inside 
two disks (of diameters independent of the wavelength and equal 
to the maximum diameter of each estimated source) are shown as a 
function of the wavelength, both for the originals and reconstructed 
sources. Again we find a good match, despite the sparse Fourier 
coverage. 

These results prove the efficiency of the adopted strategy, which 
allies redundant dictionaries, algorithmic acceleration and parallel 
implementation. This opens the possibility of providing accurate 3D 
restoration for large scale problems, involving both high spatial and 
spectral resolutions. 

The authors thank M. Vannier, R. Petrov and F. Millour for fruit¬ 
ful discussion about the use of the differential phase and R. Flamary 
for the ADMM implementation. 
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